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We study the spatial correlations of quantum fluctuations that can be observed in multi-mode 
spontaneous parametric down-conversion in the regime of high gain. A stochastic model has been 
solved numerically to obtain quantitative results beyond the stationary plane-wave pump approxima- 
tion. The pulsed shape of the pump beam and other features of the system, such as spatial walk-off 
and diffraction are taken into account. Their effect on the spatial quantum correlations predicted 
by the plane-wave pump theory is investigated, both for near field and far field measurements, in a 
type I and in a type II phase-matching configuration. 
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I. INTRODUCTION 



The spatial aspects of quantum optical fluctuations have been the object of several studies in the past [1-3]. In 
general they show up in nonlinear optical processes, typically wave-mixing phenomena which involve a large number 
of spatial modes of the electromagnetic field. Recently there has been a renewal of attention because of new potential 
applications which exploit the quantum properties of the field for image processing or multi-channel operations. 
Examples are quantum holography [4] , the quantum teleportation of optical images [5] , and the measurement of small 
displacements beyond the Rayleigh limit [6]. An overview of this relatively new branch of quantum optics, for which 
\q ' the name quantum imaging was coined, can be found in [7]. 

The process of frequency down-conversion is particularly suitable for this kind of applications because of its large 
emission bandwidth in the spatial frequency domain. We consider spontaneous parametric down-conversion (SPDC) 
taking place in a crystal with a second-order nonlinearity set in a travelling-wave configuration. In this process 
■ the photons of a high intensity pump field are split into pairs of photons of lower energy and momentum through 
O |' the nonlinear interaction with the medium. Since no signal field is injected, down-conversion is initiated only by 
, 1 „ . vacuum fluctuations that equally cover all spatial and temporal frequencies. The spontaneous fluorescence pattern 
that arises has therefore the angular spectrum determined by phase-matching conditions, and depends only on the 
&3 , linear dispersion properties of the nonlinear material. 

Recently, we used multi-mode theory in order to demonstrate that SPDC is able to display spatial quantum corre- 
lation effects in the far field zone, where the particle-like character of the generated field has its clearest manifestation 
[8,9]. In particular, assuming that the pump field is a stationary plane- wave, the theory predicts noise reduction well 
below the shot noise level for the difference in the number of photons measured from two detection areas R\ and R2 
corresponding to couples of phase-conjugate (signal and idler) modes. In other terms, the photon number measured 
over the two detection areas are identical even at the quantum level. This phenomenon finds its explanation in the 
conservation of the photon transverse momentum which is fulfilled in each elementary down-conversion process: for 
each photon detected in say area the detection of its twin in R2 is ensured by this law. 

It is important to stress that these result holds both in a low and in a high gain regime. However, only in the first 
case single photon pairs can be resolved in time by the detectors and information on spatial quantum correlations can 
be obtained directly from coincidence measurements. In this paper we shall refer more explicitly to the second case, 
where a large number of photons are emitted in each mode and their detection gives rise to continuous photocurrents. 
The observation of photon number correlation phenomena in such a regime is the aim of an experiment presently 
performed at the University of Insubria at Como. In this experiment photodetection is performed by means of a high 
quantum efficiency charged coupled device (CCD) camera, which is able to resolve photon number fluctuations that 
are below the standard quantum limit [10]. The pump field is a high power picosecond laser pulse that provides energy 
for a large number of down-converted photons, in a configuration such that the plane-wave and continuous-wave (cw) 
pump approximations are very raw. 

In this paper we present a realistic description of the system, based on a numerical model that includes the finite 
frequency bandwidth of the pump, both in spatial and temporal frequency domains. Other features of the system 
that are relevant from an experimental point of view, such as spatial and temporal walk-off, different kinds of linear 
dispersion and phase-matching (type I and type II crystals), are included in the model. It is important to investigate 
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how they affect the spatial quantum correlation phenomena predicted by the plane-wave pump theory also in order 
to identify the best conditions under which they can be observed in the experiment. A numerical evaluation of the 
far field photon number correlation function is presented in [11] in the case of a type I crystal at degeneracy. In [11], 
which treat spontaneous down-conversion within a classical framework, the shape of the pump pulse included into the 
numerical model is taken from experimental data and the obtained signal-idler correlation peak displayed between 
symmetrical point reproduce well the correlation measured experimentally 

We shall also focus our attention on the spatial correlation property of SPDC in the near field, where the signal 
and idler beams are found to exhibit quantum correlated photon number fluctuations when measured from detection 
areas that image the same portion of the beam cross-section. Twin photons are indeed generated simultaneously and 
they remain localized in a limited region of space as a long as they are observed close to the crystal. This "position 
entanglement" of the generated photon pairs can be seen as the near field counterpart of the momentum entanglement 
which can be observed in the far field. However, we shall see that for a realistic crystal length the measurement of 
near field correlation is strongly affected by propagation effects, in particular diffraction and spatial walk-off. We shall 
propose a procedure to overcome at least partially this problem. 

The paper is organized as follows. In Sec. II we briefly introduce the theoretical model used to describe SPDC 
within a classical framework. The quantum description of the system is illustrated in Sec. Ill, where a fully analytical 
treatment is developed in the framework of the plane- wave and cw pump approximation (PWPA). It is based on a 
multi-mode input-output formalism, first introduced in [12] for a type I crystals at degeneracy, which is here extended 
to a type II phase-matching configuration. 

In Sec. IV we give a qualitative description of the phase-matching mechanism that determines both the photon 
number distribution and the characteristic bandwidths of the down-converted field, illustrating thereby the differences 
between type I and type II phase-matching. 

In Sec.V we define the quantities that can be measured experimentally and that put in evidence the quantum nature 
of the spatial correlations in which we are interested. Their analytical expressions are derived within the PWPA, 
which will be used to interpret the results of the numerical model. 

The last part of the paper (Sec. VI) is devoted to present the numerical results obtained for two particular crystals 
with different phase-matching (type I and type II). The amount of correlations that may be achieved is evaluated as 
a function of different parameters that can be varied experimentally, such as the size of the pump beam waist and 
the size of the detectors. 



II. CLASSICAL DESCRIPTION OF THE PROCESS 

We decompose the electric field in the superposition of three quasi-monochromatic wavepackets (denoted with 
i?o, Ei and E 2 ) of central frequencies u>o, oj\ and ui 2 , corresponding to the pump, the signal and the idler fields, 
respectively. These frequencies are taken to satisfy the energy conservation condition w\ + uj 2 = t^o- Assuming the 
mean direction of propagation is the z direction, and denoting with x — {x,y) the coordinate vector in the transverse 
plane, we can write 

£j(2 1 ? ) t)«^(«,f 1 i) e^ z -^* + c.c. (j =0,1,2). (1) 

where kj = rijLUj/c is the wave number of wave j at the carrier frequency along the z-axis (for an extraordinary wave 
the refraction index rij depends on the propagation direction, a property leading to spatial walk-off). To simplify the 
notation we have ignored the vectorial character of the three fields, their polarization being determined by the kind 
of phase-matching condition that are met inside the crystal. 

Within the paraxial and slowly varying envelope approximation, the propagation equations for the signal and idler 
(S/I) field envelopes and the pump field envelope can be written in the form [13] 

aA Ate- iA °* (j, I = 1, 2; j ± I) , (2a) 

-aA 1 A 2 e^ az (2b) 

The driving terms on the r.h.s. describe the wave-mixing process due to the second-order nonlincarity of the medium, 
the coupling constant a being proportional to the effective second-order susceptibility \ e ff characterizing the down- 
conversion process. A = fci + k 2 — fco is the collinear phase-mismatch of the central frequency components. 
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Linear propagation is described by the l.h.s. of these equations: the terms proportional to fc' = ( ^7) and 

K' = ( ) ^ eac ^ to temporal walk-off between the different waves and group velocity dispersion respectively, 

while the terms containing the first and second order derivatives in the transverse coordinates (x, y) are responsible of 
spatial walk-off and diffraction respectively, pj indicates the walk-off angle of wave j, determined by the anisotropy 
of the crystal (the walk-off direction is taken along the y-axis). Linear losses are neglected, so that the three wave 
exchange energy but their total energy is conserved. 

In Sec. VI A we shall also consider the special case of a type I phase-matched crystal where the signal and the idler 
fields have the same polarization and are observed close to the degenerate frequency lo\ = UJ2 = ojo/2. Under these 
conditions the signal and idler fields are no more distinguishable and the down-converted field must be described by 
a single slowly-varying envelope A(z,x, t) satisfying the following propagation equation 

t + k 'f + - - = ° A ° A "^" ■ < 3 » 

which is readily obtained from Eqs.(2a) by dropping the S/I indexes j, I, which denote different polarizations and/or 
carrier frequencies in the non-degenerate case. 

In a single-pass configuration with crystal length on the order of a few millimeters, the amplitudes of the down- 
converted field remain small with respect to the pump amplitude and the nonlinear driving term in the r.h.s of 
Eq.(2b) can be neglected. The pump depletion due to down-conversion and absorption is indeed of small entity, 
unless extremely high intensity laser sources are used. We shall therefore work within the parametric approximation, 
that treats the pump as a known classical field which propagates linearly inside the crystal, while the down-converted 
fields are quantized according to the rule that are briefly illustrated in the next section. 



III. QUANTUM DESCRIPTION IN THE PARAMETRIC APPROXIMATION 

We need now to substitute the classical signal and idler fields with operators. Making the formal substitution for 
the field envelopes Aj(z, x, t) — > a,j(z, x, t) , (j = 1, 2), we impose the following commutation rules at equal z [2] 

ai {z,x,t),a){z,x',t')\ = 5 ij S{x-x ')5{t-t') , (4) 
[a,i(z,x, t),a,j(z, x" ,t')] =0 = 1,2), 

valid within the framework of the paraxial and quasi-monochromatic approximations. With this definition 

Ij(z,x,t) = a](z,x,t)a j (z,x,t) (j = 1,2), (5) 

is the photon flux density operator associated to wave j: its expectation value gives the mean number of photons 
crossing a region of unit area in the transverse plane. In the linear regime the field operators obey the same equations 
as the corresponding classical quantities. To our purposes, it is useful to introduce the Fourier transforms of the field 
envelopes with respect to time and to the transverse plane coordinates: 

aj (z,q,Q) = J^J ^=a,(z,x,t)e-^ nt (j = 1,2) . (6) 

A similar definition holds also for the Fourier component Ao(z,q, il) of the classical pump field envelope. The 
propagation equations (2a) take then the form 

d aj (zM =iSMQ)a . Mil ) (7) 



02 J^J ^=A (z,q-q',n-n')a\(z,-q',-n') (j, I = 1, 2; j ? I) 



where 



6 3 {q,n) = k' ] n + h i ;n 2 + p 3qy -^-{ q l + q l) , (3 = 1,2), (8) 
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is the quadratic expansion of kj z (cjj + Q,q) — kj around q = 0, fi = 0, and kj z (ujj + fl,q) = ^Jk 2 (ujj + Sl,q) — q 2 
denotes the ^-component of the fe-vector associated to the (q, S})j plane-wave mode. In particular the walk-off angle 
Pj can be identified as ^ calculated for q = 0, SI = 0. 

Eqs.(7) contain the convolution integral in Fourier space of the S/I field envelope with the pump field envelope. 
Within the undepleted pump approximation, the latter can be expressed as 



A a (z, q, SI) = e tSo ^ z A (z = 0, q, SI) , 

s (q, si) = k' n + h^si 2 + Poqy - ^-y x + q 2 y ) , (j = 1, 2), 



(9a) 
(9b) 



the z = plane being taken at the input face of the crystal. In the following we shall assume that the pump pulse 
has a Gaussian profile both in space and time, of beam waist wq and time duration r at z = 0: 



In Fourier space we have then the expression 

A (z = 0,q,Sl) = 2V2-^ e -^M/^V^o 2 



where 



SqfiSuJo 



Sq Q = 2/w , Sluq = 2/t 



(10) 

(11) 
(12) 



denote the bandwidths of the pump in the spatial frequency domain and in the temporal frequency domain respectively. 
Let us now consider the limit of the PWPA approximation, in which wo and To tend to infinity and 

A (z, q, SI) - (2n) 3 / 2 A p 6(q)6(Sl) . (13) 

Under this condition Eqs. (7) couple only pairs of phase-conjugated modes (q, Sl)i and (—q, — S}) 2 and can be solved 
analytically. The unitary input-output transformations relating the field operators at the output face of the crystal 
a° ut (q, SI) = a,j(z — l c , q, SI) to those at the input face a*"(<7, SI) = aj(z = 0, q, fl) take the following form 

a° ut (q, SI) = Ui% Sl)a\ n (q, SI) + V^q, S})a 2 n \-q, -SI) 
a° ut (q, SI) = U 2 (q, Sl)a 2 n (q, SI) + V 2 (q, Sl)a? \-q, -SI) 

with 

A(q,Sl) . 



(14) 



U\{q, SI) = cxp 
Vi(q,Sl.) = cxp 
U 2 (q, SI) = cxp 
V 2 (q,Sl) = cxp 



.5 1 (q,n)-6 2 {-q,-Sl)-A Q 
i L 



cosh(r(g, si)i c ) 



'2T(q,Sl) 



sinh(r(g, Sl)l c ) 



(15) 



. 6i{q, SI) - 6 2 {-q, -SI) - A , 



. S 2 (q,^)-5 1 (-g,-Sl)-A l 
i -z W 



.S 2 ( q ,Sl)-5 1 (-q 7 -a)-Ao J 

% _ If 



A(-q, -SI) 

cosh(r(-g, -Sl)l c ) + » 2r( _ :_ 0) sinh(r(-q, -Sl)l c ) 
sinh(r(-£ -Sl)l c ) , 



T(-q, -SI) 



and 



A(q,Sl) = A + <5i(9,fi) + S 2 (-q, -SI) w k lz (q,Sl) + k 2z (-q,-Sl) - k , 
cr p = oA v . 



(16a) 

(16b) 
(16c) 



It is important to note that the gain functions Uj and Vj given by Eq. (15) satisfy the following unitarity conditions: 



|^(g,ft)| 2 -|W^)l 5 



(j = l,2) 



u^zfyVii-z-n) = u^-^-s^v^si) , 

which guarantee the conservation of the free-field commutation relations (4) after propagation. 



(17a) 
(17b) 
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IV. MEAN INTENSITY DISTRIBUTION 



In the following we shall consider measurements either in the near field or in the far field zones of the nonlinear 
crystal. In order to simplify the notation we shall omit the explicit dependence of the fields on the z coordinate: when 
specification is explicitly needed, the measured quantities will be labelled with 7r or tv', which will denote the near 
field and the far field detection planes respectively (see scheme of Fig. la). The analytical results given here and in 
the next sections are all obtained within the PWPA; on the one hand they generalize those illustrated in [9] for a type 
I crystal at degeneracy to a type II phase-matching configuration, on the other hand they provide a good starting 
point to interpret the results of the numerical model that includes the pulse shape and the finite cross section of of 
the pump beam. A more general input-output formalism which goes beyond the PWPA is developed in appendix A. 

With a stationary and plane-wave pump the near field intensity distribution in the output plane of the crystal 
clearly does not depend on x and t, because of the system invariance with respect to translation in time and in the 
transverse plane. Using input-output relations (14) and recalling that the input fields is are the vacuum states, we 
obtain easily 



(Ij(x, t)) 7 



The function \Vj(q, f2)| 2 gives the contribution of mode (q, to the total photon flux of beam j, and is usually 
referred to as its spectral gain. On the other hand, in the far field plane n' the spatial Fourier modes are resolved 
spatially and the photon distribution reflects the g-dependence of these spectral functions. From the expression of 
T(q, Q) given by Eqs. (16), we see that down-conversion occurs most efficiently for the modes satisfying the condition 
A(<f, Q) < 2<r p . Using Eqs. (16b) and (8), the phase-mismatch accumulated during propagation can be written in the 
form 

A(q, Si)l c = A l c + sign[fci + p 2 q y S^±Sv (19) 

where we assumed that the signal wave is ordinarily polarized, so that p\ = 0, and we introduced the parameters 



\J lc \k'i k' 2 \l c y {k\ -\- k 2 )l c 

where k = 1k\kil(k\ + fe). They determine the characteristic bandwidths of SPDC both in the temporal frequency 
domain and in spatial frequency domain. In the type I phase-matching configuration we will consider in Sec. VI A, both 
the signal and the idler waves arc ordinarily polarized and are observed close to degeneracy, i.e. for ui\ = ui 2 = <^o/2- 
In this special case the temporal bandwidth is determined by fio = f2 ' = vVI^iMj smce ki(ui) = k2(u>) implies that 
Q' Q = co. Far from frequency degeneracy the emission spectrum has a much narrower bandwidth, on the order of 
f2o = which is about 2 to 3 orders of magnitude smaller than f2 ' (for a typical crystal length of few millimeters). 
On the other hand, in type II crystals the signal and idler waves are characterized by different polarizations and 
frequency dispersion relations, so that Slo = remains finite even for uj\ = oj 2 . The spatial bandwidth q gives the 
range of transverse wave- vectors for which the gain spectrum \Vj(q, f2)| 2 is close to its maximum value, sinh 2 [a p l c ] , 
when the j th field is observed at a given frequency ujj + £1 (the gray region shown schematically in Figs.l b and c for 
fl = 0). Wc remark that 1/ ' p2l c and qo are about the same order of magnitude as long as l c remains in the millimeter 
range. 

For definiteness we assume that the far field is observed in the focal plane of a thin lens of focal length / which 
performs the Fourier transformation of the field from the output face of the crystal (the so-called / — / system). The 
field operators in the focal plane tt' at z = l c + 2/ (see Fig. la), which we denote with 61,2(2?, £), are related to those 
in the output plane of the crystal by the following Fresnel transformation 

bj(x,t) = I dx 'hji^x^af 1 *^ ' ,t) , (21a) 

hj(S,S') = fi e -^ ' , (j = 1,2) , (21b) 

where Xj = 2izc/ujj 1 (j=l,2), are the free-space wavelengths corresponding to the carrier frequencies. Using the input- 
output relations (14) and unitarity relations (17) we can evaluate the mean intensity distribution of the two fields 
with the following approximate expression 
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« -(Jr / fV,(^)| 2 , (i = 1,2) . (22) 

where we introduced the barred gain functions defined in real space 

Uj(x, SI) = U, (p-x, n) , Vjix, SI) - V- f , Si) , (j = 1, 2), (23) 



Xjf" J' JV ' ' ' W 

and = (Xjf) 2 /S A , (j=l,2), denotes the resolut ion areas in the far field plane at the S/I wavelengths, with Sa 

being the area characterizing the dimension of the system in the transverse plane. As it is shown in [9], Eq.(22) 
can be obtained by assuming that a pupil of area Sa 3> i s P ut 011 the crystal exit face. Assuming that the 

transverse dimensions of the crystal are large compared to the pump waist, Sa can be identified with the effective 
cross section area of the pump beam. Eq.(22) represents a good approximation provided the pump beam shape 
changes negligibly during propagation in the crystal and behaves therefore as a plane-wave. This happens when the 
Rayleigh length characterizing the Gaussian pump beam divergence, = itwq/Xq, and its analogue characterizing 
dispersion, z® isp = Tq/2/cq, are much longer than the crystal length l c . The same conditions can also be written in 
terms of the pump spatial and temporal bandwidths (12) as 

(24) 

(jo "o 

At the considered carrier frequencies oj\ and L02 the gain functions Vj (q, SI = 0) are maximal and perfect phase- 
matching is achieved when the equations A(±<f, SI = 0) = are satisfied, with the plus sign for field 1, and the minus 
sign for field 2 (see Eqs.(15) and (16b)). More explicitly, as can be seen using expression (19), they can be written as 

^|+^±Ip 2 Z c ^ = A Z C + Qp2/ C «b) • (25) 

Provided that A Z C > —jpfycloi we have therefore two circles of radius qa and centered at (q x = 0, q y = ±qc), with 

Qc = ^Pihql = ^kp 2 , (26a) 



q« = ^.\/a /, • % \ /•A,, • - { (2(il») 



They are plotted in Fig. lb, c, respectively for a type II and a type I phase-matching configuration. Modes close to 
these circles within qo (gray annuli in the figure) give a non negligible contribution to the down-converted field. In 
the detection plane they give rise to characteristic couples of rings which have been observed in many experiments on 
spontaneous parametric down-conversion (see e.g. [15-17]). It should be stressed that, without any spectral filtering, 
emission occurs on a very wide range of wavelengths and emission angles, as allowed by the phase-matching conditions 
(see e.g. [16]). However, from an experimental point of view, a particular couple of rings can always be selected with 
the use of frequency filters centered at the chosen frequencies U\ and uo 2 = uq — u)\- Noting that there is the following 
mapping between the spatial frequency plane and the far field plane: 

Xif 

(qx,q y ) -» -^r(qx,q v ) for field 1, (27) 
(qx,q y ) -> ^-(qx,q y ) for field 2, 

it is easily seen that the ring radii, x^ 2 " 1 = A *' 2 ^ qn, and their distance from the z-axis, y^ 1 '^ = Xl 2 ^ J ' qd are generally 
different except when observation is performed at frequency degeneracy (i.e. for Ai = A 2 ). In a type I phase-matching 
configuration these rings are concentric, since there is no spatial walk-off between the two fields (p 2 = 0) and the 
radial symmetry of the system is preserved. 

Fig. 2 and 3 illustrate the kind of far field patterns that can be obtained at frequency degeneracy in a type I and a 
type II crystal, respectively, assuming the pump field has a large beam waist (in practice, the condition 5qo <C qo must 
be fulfilled). They are obtained by numerical integration of the classical field equations (2), with a white input noise 
which simulates the vacuum fluctuations that trigger the process, as it will be described in the next sections. These 
2-D simulations do not include the temporal dimension and cannot be used to obtain quantitative results; however, 
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they provide some insight on the spatial features of the far field patterns generated in single pump shots, as they 
could be observed experimentally by using a narrow band filter at £1 = 0. 

In the type I phase-matching case, being at frequency degeneracy, the two rings merge into one that contains both 
the signal and idler modes. It should be noted that qn and qc do not depend on the crystal length whereas qo scales 
as lj\JT c , so that the thickness of the rings at a given temporal frequency is larger for a shorter crystal. 

V. NEAR AND FAR FIELD CORRELATIONS 

We now define explicitly the quantities that can be measured in an experiment in order to put in evidence the 
S/I correlations in the spatial domain we are investigating. We assume that the signal and idler beams are spatially 
separated in the detection plane and are measured over two detection areas which we denote with R\ and R 2 . 

In the far field, correlations find their origin in the conservation of the transverse momentum of the generated 
photon pairs. Therefore, in order to find maximal correlation, R\ and R 2 must correspond couples of phase-conjugate 
modes, such as those indicated with the black squares in Fig. 1 b,c. For simplicity, in order to avoid the heavy 
notations which arise if Ai =/= X2, we shall restrict our analysis to the frequency degenerate case, indicating with A 
both Ai and A 2 . Phase-conjugate modes are then mapped by the lens into symmetrical points in plane tt' according 
to relation (27) and Ri and R 2 must be taken symmetrical. 

On the other side, near field correlations arising from the position entanglement of the twin photons are expected 
to be observed if i?i and R2 occupy the same region of the near field plane. In practice a type II phase-matching 
configuration should be considered, so that the use of a polarizing beam-splitter and lens systems allows the imaging 
of the S/I near fields on two physically separated detection planes (see detection scheme illustrated in Fig. 4). 

If the detectors are the pixel of a CCD camera, as in the experiment described in [16], they do not allow any 
spectral measurement due to the very low resolution power of the device in the time domain. They simply measure 
the total number of incoming photons down-converted in each single pump shot and the measurement time Td can 
be identified with the pump pulse duration. We introduce therefore the operators corresponding to the number of 
photons collected by the two detectors in the finite time window [—Td/2, Td/2] 

Nj=[ dx I 1 dtlj(x,t) (i = l,2), (28) 

JRj J-T d /2 

The measurable quantity which is capable of displaying the quantum nature of the photon number statistics in the 
spatial domain is the variance of the photon number difference, 7V_ = N\ — N 2 , which can be written in the form 

((5N_) 2 ) = (N+) + (: (8N,) 2 :) + (: (SN 2 f :) - 2(5^8^) . (29) 

SNj = Nj — (Nj) and 5N- = N- — (N_) denote the photon number fluctuation operators associated to Nj and iV_ 
and the colon ":" denotes normal ordering (n.o.) for the expectation values. In Eq.(29) the shot noise contribution, 
i.e. the total number of photons intercepted by the two detector (N + ) = (Ni) + (N 2 ), has been explicitly separated 
from the term which describes the field correlations. We define 

T /2 T I '2 

(: SNiSNj :) = / dx [ dx ' [ dt [ dt'G t3 {x,t,x ' ,t') , = 1,2), (30) 

JRi JRj J-T d /2 J-T d /2 

where 

Gij(x,t,x',t') = (: Ii(x,t)Ij(x',t') :) - (h(x, t))(Ij(x ', t 1 )) , = 1,2). (31) 

are the n.o. self- and cross-photon number correlation functions of the S/I beams. Notice that in the non-degenerate 
case (type II or type I far from frequency degeneracy) we have (: I\(x, t)I 2 (x ' , t') :) = (h(x, t)I 2 (x 1 , t')). 

We now focus on the analytical results that can be deduced from the PWPA. We shall consider explicitly only the 
case of type II phase-matching. The case of type I, at least in the far field, can be described with a similar treatment 
and has been already discussed in [9] . Assuming that the detection time Td is large compared to the coherence time 
T co h = fig ' as ^ 1S usually the case, we have 

(: SNiSNj :) w T d [ dx j dx' G l3 {x, x' , O = 0) = 1,2) . (32) 

JR t JRj 

where Gij (x, x ', Cl = 0) is the Fourier transform of the function (31) with respect to t — t' (notice that in a cw regime 
this function depends only on t — t'). Next, the Gaussian character of the field statistics allows to express fourth order 
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correlations in terms of second-order correlations (see App.A, Eq. (A5b)) In this way the photon number correlation 
defined by Eq. (30) in a plane z can be written as the photon number correlations defined by Eq.(30) in a plane z 
can then be written as 

{:&Ni&N r ) t =T d I dx I dx" j ^\Y^\x,x',ir)\ 2 , (i,i = l,2), (33) 



where 



T ( j f(x,x',ri) = JdTe- inT (a](z,x,t + T)a J (z,x',t)), (j = 1,2) (34a) 

r[ z 2 \x,x ',Q) = J dTe- mT ( ai (z 7 x,t + T)a 2 (z 7 x ',*)) , (34b) 

are the only correlation spectra of the S/I fields which do not vanish when the input field is in the vacuum state. We 
assumed here implicitly that propagation in free space occurs without losses. 
In the far field plane n' the self and cross correlation functions are given by 

r^p(x,x',n)=5(x-x')\V j (x,ft)\ 2 , (.7 = 1,2), (35a) 
r^\s,S' ,0) = -6(x + x ')Ui(x,n)V 2 (-x,-n) . (35b) 

as it can be inferred using the Fresnel transformations (21). In this case, both correlation functions display a delta-like 
peak, located at x ' = x for the self correlation, and at x ' = — x for the cross correlation. The delta-like character of 
the correlations derives from the unphysical assumption that the transverse dimensions of the system are infinite. In 
[9] which deals only with far field correlations the transverse size of the system were taken into account a posteriori 
by considering a finite aperture Sa placed at the crystal output face, assuming a condition equivalent to (24) is met. 
With this approach, we found that far field correlations are localized within the resolution area determined by the 
size of this aperture (that is Sdiff — (A/) 2 /Sa with the / — / lens system). This procedure also eliminates the 
cumbersome divergencies arising from the singularity of the spatial delta-functions, allowing the formal substitution 
S(x = 0) — > 1/ Sdiff as has been done in Eq.(22). By integrating over two symmetric detection pixels of area much 
larger than this resolution length we obtain 

(: (SN,) 2 :>„, = (: (5N 2 f :)„, = J— f f dx Sl)| 4 , (36a) 

Sdiff J 2?r J Rl 

(5N 1 6N 2 ) 7V ' = I ^ I dx \U 1 (x,n)V 2 (-x,-fl)\ 2 . (36b) 

Sdiff J 2tt J Rl 

In order to obtain these expressions, both unitarity relations (17a) and (17b) must be used, together with the fact 
that the integration area R\ and R 2 are taken symmetric with respect to the origin. We then easily get 

(: OWx) 2 + (: (^ 2 ) 2 - 2<<W 1 <W 2 ) 7r , = -_L- f ^ f dx Fi(x,0)| 2 , (37a) 

Sdiff J 2vr J Rl 

= -(N+) w> , (37b) 

which implies (((5A_) 2 ) 7r / = 0, as it follows from Eq. (29). A more rigorous approach to the issue concerning the finite 
resolution of the system is given in appendix B, where we derive an approximate solution of the propagation equations 
which include the finite pump dimensions in the limit (24). It is shown that the width of the far field correlation 
peaks are indeed on the order of x d iff — (\f /2ir)5qo, the resolution length determined by the pump beam waist Wq. 
In the near field plane tt (at z = l c ), using Eqs.(14) and (4) we obtain the expressions 

r%\z,Z',0) = J ^e-^-Oiv^fi)^ (j = 1,2), (38a) 

T[l\x, x',n) = J J^ e *(*-* ')[/!(£ n)V 2 (-q, -SI) , (38b) 

which depend only on the offset between the two points, x — x ', as a consequence of the invariance of the system 
with respect to translations in the transverse plane which follows from the PWPA. Provided that the typical scale of 
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variation of the function appearing under the integrals in Eqs. (38a) and (38b) is q 0} we expect that these correlations 
are localized in a region of size 



x co h = l/9o ~ \[hjk , (39) 

a quantity which can be identified with the transverse coherence length of the down-converted fields. This finite 
correlation length comes from the spread out of the generated photons due to diffraction, which increases proportionally 
to the square root of the propagation distance; we can therefore expect that the detection areas must be larger than 
this coherence area in order to measure good correlations in the near field. 

We incidentally note that the cross correlation function Ti 2 displays a localized peak, that in the near field is located 
at x ' = x (Eq.(38b)), while in the far field is located at x ' = —x (Eq.(35b)). The x <-> x ' correlation in the near 
field reflects the entanglement in position of the twin photons, while the x <-> —x 1 correlation in the far field comes 
from their entanglement in momentum. 

For what concerns the near field, let us consider more in general what happens in a generic plane of coordinate z 
close to the crystal output face. As it is shown schematically in Fig. 4, we consider a measurement in which the signal 
and the idler fields are separated with a polarizing beam splitter placed beyond the crystal. The two lenses L and 
L' put in the signal and idler arms performs the imaging of plane z onto two distinct detection planes (note that we 
have here two 2/ — 2/ lens systems, while in the case of the far field measurement considered previously we had a 
single / — / system). The S/I photons are collected by means of two square pixel detectors R\ and R 2 , centered at 
the positions x\ and x 2 respectively. The propagation from the crystal exit face z = l c to the detection planes can be 
described by a Fresnel transformation of the form (21a) with the kernel 

h{x,x') = ~\ e Tral*-*'l 2 (40) 
\(z — l c ) 

where inessential phase factors due to the presence of the lenses have been omitted. By using this transformation 
inside Eqs. (34a) and (34b) we can calculate explicitly each term on the r.h.s. of Eq.(29). By performing explicitly 
the integration over the square pixel areas in Eq.(33), we obtain 

(: (SN 3 ) 2 :) z =T d J dqj dq'H n (q,q') J tyflW, V\ 2 , (j = 1,2) (41a) 

(SmSN^^Td J dqj dq'H 12 (q,q>) J ^ U l (q, n)V 2 (-q, -flpl (q ', H)V* 2 (-q ', -fi) , (41b) 

(N + ) z = 2T d d*J J £L J ^\VMm 2 (41c) 

where the functions depend on the square pixel size d and on their relative positions through the relations 

Tl„ „>\Jl Yin -n'\rll 

(42a) 



~(q x - q' x )d~ 




'{Qv ~ <l'y) d ~ 


2 


sine 2 


2 



Hu{q,q') = (d/2n) 4 sine 2 

Hi2{q,q") = e- ih ^( q2 - q '^ +i tf-t'>( Sl -^H n {q,q'). (42b) 

We first note that if the detection areas are reduced well below the coherence area x 2 coh , the fluctuations of iV_ 
approach shot noise. Indeed, in the limit d <C x co h we can replace Hu(q,q') with Hu(q,q) = d i /(2ii) A in Eqs. (41a) 
and (41b), from which it can be verified that the correlation terms scale as d jx coh while the shot noise contribution 
scale as d 2 /x^ oh (to evaluate this scaling it should be noted that for a fixed f2, the area in (f-space where the gain 
functions are not negligible is on the order of q$ = ^/x 2 coh ). On the other hand, if the detection areas are large enough 
with respect to x 2 oh , the substitutions 

Hv(q,q')^J^5(q-q') (43) 

can be used for evaluating both (41a) and (41b); this leads to vanishing fluctuations in the measurement of iV_ 
as for the far field case (again, unitarity relations (17) must be used in order to obtain this result). However the 
condition " large enough" is more stringent for cross-correlation than for self-correlation. Let us first consider the case 
z = l c , and x\ = x 2 - The function appearing under the integral in (41a) is always positive. By contrast, the function 
appearing under the integral in (41b) is an oscillating function, which becomes positive only for an infinite pixel size, 
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when the limit (43) is strictly achieved. This feature tends to lower (SNiSN2) z with respect to (: (SNi) 2 :) z and 
(: (6N2) 2 :) z . As a consequence of this behavior the fluctuations of iV_ will therefore exceed shot noise, as it can be 
easily inferred from expression (29). As a matter of fact, as we shall see in Sec. VIC, x\ — X2 and z — l c are not the 
better choices to minimize the fluctuations of N- , and special care in the positioning of the detectors is necessary in 
order to compensate both the effect of diffraction and spatial walk-off between the S/I fields, which are included by 
the phase of the function U\(q, Q)V2(— q, — fi). 

Finally we note that if the losses of the detection process are taken into account, the ideal result {(5N-) 2 ) = must 
be replaced with 

((5N-) 2 ) = V (1-V) (N+) , (44) 
■q denoting the finite quantum efficiency of the detectors. 

VI. NUMERICAL RESULTS 

We now present the results obtained from the numerical model that includes the effects of the finite pump. The 
quantum averages in which we are interested (i.e. mean photon numbers and photon number correlations) are 
evaluated through a stochastic method based on the Wigner representation. With respect to other representations 
in phase space, the Wigner representation presents the advantage that the c-number stochastic equations equivalent 
to the equations for the field operators (2) do not contain Langevin noise terms (because of linearity and absence 
of dissipation) and are therefore identical to the classical propagation equations (7). The statistical character of the 
quantum fields is wholly contained in the stochastic input field. We therefore proceed as follow: 

1) We generate the input field with the appropriate phase-space probability distribution, which is a Gaussian white 
noise with zero mean, corresponding to the vacuum state in the Wigner representation [19]. 

2) We perform the numerical integration of Eqs.(7). We use a split-step algorithm [18] which integrates separately 
the terms describing linear propagation and the term describing the wave mixing process: the former are integrated 
in Fourier space, the latter in real space. 

3) The obtained output fields are used to evaluate the correlation functions of interest. The procedure must be 
reiterated a sufficiently large number of times, so that the stochastic averages performed becomes good approximations 
to the corresponding quantum expectation values. Furthermore, some corrections are usually necessary in order to 
convert them to the desired operator ordering (the Wigner representation yields quantum expectation values of 
symmetrized operator products). 

A. Far field correlation in a type I crystal at degeneracy 

We first focus our attention on far field correlations that can be observed in a type I crystal with emission close to 
degeneracy (the case illustrated in Fig. 3). Since there is no spatial and temporal walk-off between the signal and the 

idler modes, the most significant parameters in play are the spatial bandwidth q = \J k/l c and the the ring radius 

qji = v fcAo, which are determined by the crystal length and the collinear phase-mismatch parameter (see Eqs.(20) 
and (26b)). To a large extent, the analysis that we have performed docs not depend on the particular type I crystal 
that is considered. We shall however refer to the specific case of LBO (LiB 3 5 : lithium triborate) set in the type I 
configuration described e.g. in [15,16]. The pump field operates at Ao = 532 nm and propagates in the XY crystal 
plane forming an angle with the V-axis close to 11°, for which collinear phase-matching at degeneracy is achieved. In 
order to evaluate the characteristic bandwidths and walk-off parameters given by Eq. (20), we have used the Sellmeier 
dispersion relation coefficients which can be found in [13]. Considering a crystal length of 5 mm, the coherence time 
is T co h = ^0 1 ~ 0-01 ps. Assuming the pump pulse has a duration of about 1.5 ps, as in the experiment described 
in [16], the ratio Sujo/^Io is as small as 10~ 2 . The analytical results obtained within the plane-wave and cw pump 
approximation are therefore expected to provide good insight as long as the ratio Sqo/qo remains small compared 
to unity. The temporal/spatial walk-off between the signal and the pump field are negligible unless the pump pulse 
time/beam waist are exceedingly small, a situation we do not consider, since it is too far from the ideal plane wave 
pump limit (it would prevent the observation of any spatial quantum correlation effects). 

Fig. 5 displays the intensity distribution obtained numerically from a single integration of the propagation equations, 
for increasing values of the pump beam spatial bandwidth Sqo = 2/uuq, both in the near field and in the far field planes. 
As already mentioned, these simulations with two transverse dimensions do not include the temporal dimension and 
are only meant to provide some insight on the spatial features of the down-converted field generated in a single pump 
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shot. In the near field, the intensity follows the Gaussian profile of the pump and displays a noisy spot pattern with 
a characteristic wavelength n/qR. In the far field, the intensity peaks (white spots in the figures) always appear in 
symmetrical pairs as in the plane-wave pump case illustrated in Fig. 3 b. However, they become broader and broader 
as the pump beam waist wo beam is reduced. Their size in the observation plane n' is on the order of the resolution 
length imposed by the finite transverse size of the pump beam waist, i.e. x<nff — ^&Zo- The following quantitative 
evaluations of the amount of correlation include instead the temporal dimension but consider only one transverse 
dimension in space, since a full 2-D calculation (2 transverse dimensions + time) would have required an exceedingly 
long CPU calculation time for our computer. 

Some insight can be gained by looking at the n.o. photon-number correlation (: SN(x)SN(x ') between two 
pixels centered at x and x ', defined by an expression similar to Eq.(30) which refers to the type II phase- matching 
configuration. The vectors x and x ' denote here the positions of the two pixel detectors in the transverse plane 
(typically two pixels of a CCD camera), their size being determined in the simulation by the spatial step of the 
numerical grid. In the degenerate case considered here the indexes i,j are dropped since the signal and idler fields 
are not distinguishable. Integration in time is performed over an interval Td which is taken larger than the pump 
pulse time, To = 1.5ps, so that all the down-converted photons generated by a single pump pulse are collected by the 
two detectors. Its 1-D (one transverse dimension + time) numerical evaluation is plotted in Fig. 6 as a function of x, 
keeping x' fixed at x' — (A/ /2n)qR, where the gain is maximum, and for different values of the ratio Sqo/qo- The 
x-coordinate is normalized to x = ^qo, the spatial scale of the photon number distribution in the far field plane. 
According to these simulations, the widths of the two peaks is on the order of Xdiff, the resolution length imposed 
by the pump beam transverse dimensions. The correlation peak at x = —x', left in the plots of Fig. 6, is always 
higher than the correlation peak on the right at x = x'. In a similar way as it was shown in Sec.V (see Eqs.(29) and 
Eqs.(37)), this particular behavior of the n.o. correlation function indicates the possibility that the fluctuations of 
N- = Ni(x) — N2(x') vanish when measured from two symmetrical pixels (i.e. taking x = —a?'). Indeed, for two 
disconnected pixels we have 

((SN_) 2 ) = (N+) + (: (SN(x)) 2 :) + (: (SN(x')) 2 :) - 2{5N(x)SN(x')) (45) 

where (N + ) = (N(x) + N(x')) represents the shot noise for 7V_; therefore, since ((5N-) 2 ) is always a non negative 
quantity, the following inequality holds: 

2(5N(x)SN(x')) - (: (SN(x)) 2 :) - (: (SN(x')) 2 :) < (N+) . (46) 

In particular, when the equality sign holds in Eq. (46), the maximum amount of correlation between the points x and 
x' is achieved, implying thus a complete suppression of the noise of N-. 

However the amount of noise in N- depends on the actual size of the detectors used to probe the correlation. We 
have evaluated numerically the variance of N- for two symmetric detection areas, by varying the size of the detectors 
as could be obtained in practice by grouping several pixels of a CCD. The results are shown in Fig. 7, where the 
different lines correspond to different values of the ratio Sqo/qo, that is to different values of the pump beam waist. 
Fluctuations are well below shot noise when the detector size d is larger than Xdiff, that is for d/xo > 8qo/qo, which 
in turn implies that the detection size must be larger than the width of the correlation peaks. These simulations show 
that the localized character of the correlation predicted by the plane-wave pump theory is well preserved as long as 
the pump beam waist is not too small. Only in the worst case considered, with Sqo/qo — 0.5, the noise reduction 
factor ((iV_) 2 ) 7r //(A^ + ) 7I -/ is is never close to zero unless the detectors cover the whole widths of the the ring pattern. 
It should be noted that a further increase of the ratio 5qo / qo would lead to single- mode emission: the inverse of these 
ratio provides indeed an estimation of the number of spatial modes that are efficiently amplified when the field is 
observed at a fixed temporal frequency. 

B. Far field correlation in type II crystals 

For the case of a type II phase-matching configuration, we shall consider explicitly the system described in [16]: 
a 1.5ps high intensity laser pulse is injected in a 4 mm long beta barium borate (BBO) crystal cut for type II 
phase- matching. In the example we consider the pump is oriented at an angle close to 48.2° with respect to the 
crystal axis and SPDC is observed around the degenerate wavelength Ai = A2 = 704nm with a lOnm interference 
filter. For the chosen parameters, the radii of the two rings q^ vanish, which is the situation illustrated in Fig. 2c. As 
for the type I configuration, we investigated on the momentum correlation that can be observed in the far field plane 
7r', by considering two symmetrical detection areas. The variance of A_ normalized to shot noise is plotted in Fig. 8 
as a function of the detector size d and for different values of the pump beam waist. The result is similar to the one 
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obtained for the type I crystal configuration illustrated in Fig. 7: fluctuations are well below shot noise only if d is 
larger than the characteristic resolution length of the system = ^^q . We note that the number of temporal 

modes that are amplified by the crystal is much lower than in the degenerate type I configuration, since in this case 
temporal walk-off between the signal and the idler fields reduces drastically the emission bandwidth Q - Indeed, 
for the picosecond pump pulse we considered here Sluq/^o is in the order of the unity for the type II configuration, 
against the 0.02 value found for the type I LBO crystal at degeneracy, for which temporal walk-off is not present. 
The numerical simulations show that this feature of type II phase-matching does not affect spatial correlations, but 
simply lowers the number of generated photon pairs per pulse. 



C. Near field correlation in type II crystals 

A main advantage of the type II configuration lies in the fact that the signal and the idler fields have different polar- 
izations and can therefore be manipulated more easily. In particular, it is possible to measure their mutual correlation 
in the near field after they have been physically separated by a polarizing beam splitter, as shown schematically in 
Fig. 4. The lenses L and L' shown in the figure simply perform the 2/ — 2/ imaging of the "near field plane" 7r onto 
the two detection planes. For the moment, we only assume that plane 7r is located at some coordinate z inside the 
crystal. The S/I field self- and cross-correlation functions, which in the plane- wave pump limit have the expressions 
(38a) and (38b), display pronounced peaks for x ' — x. In particular, the cross-correlation peak of r 12 (x, x ',51) 
describes the position entanglement of the S/I photons, which are generated in pairs in the same region of the crystal. 
The width of the peaks is on the order of the coherence length x co h = l/<Zo defined in Eq.(39). It reflects the spread 
out of the generated photons due to diffraction, which increases proportionally to the square root of the propagation 
distance. In addition, depending on the phase-matching conditions, twin photons can be emitted non collinearly with 
an aperture angle on the order of a_R = 2qn/k 1 qn being the radius of the rings in Fourier space given by Eq.(26). This 
introduce a further indeterminacy on the order of o.rI c /2 = qR/qox co h in the relative positions of the twin photons 
measured in the near field. We expect therefore that the two macroscopic fields display identical fluctuations when 
observed from the same region of the near field plane, provided that the detection areas i?i and R 2 are larger than the 
indeterminacy introduced by these propagation effects. However, for crystal lengths on the order of a few millimeters 
this indeterminacy can be as large as several tens to hundreds of microns, so that a substantial portion of the two 
beams must be intercepted in order to measure significant correlation in the quantum domain. We shall see with a 
specific example that this difficulty can be at least partially overcome. 

The role of spatial walk-off is more subtle to determine. As pointed out in [24], the Poynting vectors of the phase- 
matched modes, which determine the photon fluxes of the signal and the idler beams, generate two cones which have 
the same axis inside the crystal. To verify this point explicitly we consider the Poynting vectors associated to two 
particular modes of the signal and idler fields at ft = with transverse wave- vectors qi and q 2 . Within the paraxial 
approximation their directions are determined by the (two dimensional) angles 

Qi - 9*2 . „ 

«! = -, ^ = V2 +P2, (47) 

where p 2 = (0, — p 2 ) indicates the walk-off direction of the idler field (the minus sign is due to the fact that in the 
chosen reference frame the walk-off is oriented opposite to the y-axis, po and p 2 being assumed to be positive) . For the 
modes propagating along the axis of the signal and idler cones, we have q\ = (0, — qc) and q 2 = (0, qc) with qc — \kp 2 , 
from which we see that the corresponding Poynting vectors are collinear with 6l\= 6t 2 = (0, — k ^+ k Pi)- As a result, 
if the near field is measured directly on the crystal output face z — l c , walk-off do not contribute to tire indeterminacy 
in the relative position of the twin photons. It should be stressed, however, that free propagation beyond the crystal 
occurs at angles which are simply proportional to the transverse wave-vector of the phase-matched modes, since in 
free space the Poynting vector and the fc-vector directions again coincide. Outside the crystal the signal and idler 
emission cones are therefore oriented along different directions with an aperture angle ac = 2qn/k^ = -g^jp 2 , as 
illustrated schematically in Fig. la, and give rise to separate rings in the far field (k^ — 2n/Xi — 2n/\ 2 denote here 
the wave number of the S/I fields in free space at the career frequencies). In the simulation illustrated in Fig. 9, we 
consider the type II BBO crystal in the same conditions described in Sec. VI B. The near field coherence length is 
x C oh — 16.6^m and q^ — 0. The plot displays the noise reduction factor ((iV_) 2 ) 7r /(iV + ) 7r , evaluated numerically as a 
function of the 1-D detector size d. If the near field observation plane 7r coincides with the output face of the crystal at 
z = l c (white circle), we see that the fluctuations are significantly reduced only when d is about 15 times larger than 
Xcoh. The improved result (black squares) has been obtained by imaging onto the detection planes a plane inside the 
crystal at z — l c — Az, rather than the crystal output face. Furthermore, the array of pixel detectors in the signal and 
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idler arms are shifted with respect to each other by a distance Ay in the transverse direction of walk-off. Notice that 
this is a 1-D simulation, and the quantities Ay and Az correspond to the quantities x\ — x 2 and l c — z, respectively, 
which appears in the definition (42b)) of the function i?i2 (<?*, q' )• Az and Ay are chosen in order to minimize the 
dependence on q and q 1 of the phase of the integrand in Eq.(41b), maximizing in this way the photon- number cross 
correlation (SN1SN2). This is achieved by taking 

tanhcrp/ c m + n 2 . 

Az op t = f l c , (48a) 

2a p l c 2nin 2 

tanh a v l c , r ^ 1 

Ay op t = - , P P2h ■ (48b) 

Indeed, with these values of Ay and Az the phase factor of the function Hi2(q, q') defined by Eq.(42b) nearly cancels 
the phase of the gain function product appearing in the r.h.s. of Eq.(41b), as can be verified by using the approximate 
expression 

aigpi&Sim-Z-Witf'MVti-t'i-to)] « Y C [A(q,n)-A(q\n)]l c , (49a) 



2aJ 



„2 „/2 



tanhcr p Z c t q" — q' 



P2{q" v ~ €y) 



9o 



(49b) 



which holds in the high gain region of the spatial frequency plane. These shifts of the detection S/I planes are 
necessary in order to minimize both the effects of diffraction and of spatial walk-off. 

The plot of Fig. 10 displays the variance of A_ normalized to shot noise in the (Az,Ay) plane, as calculated 
from Eqs.(41). In this example, d is only twice the coherence length, x co h = 16.6/im, while Ay op t — 47.5/im and 
Az opt — 407/im (with a p l c = 3). The fluctuations are well above shot noise everywhere except in the narrow diagonal 
region around the point (Az opt , Ay opt ), where ((SN-) 2 ) 7T /(N + ) 7T = 0.3. Outside this region, the self-correlation 
become much larger than the cross-correlation and the variance of A_ rapidly exceeds the shot noise level. This 
means that for such small detectors a highly precise imaging of the near field determined by Eqs.(48) is therefore 
necessary in order to observe some noise reduction effect. 

Although the previous results rely on a detailed description of the relative phases of the signal and idler fields, we 
can give them a more intuitive explanation based on the particle picture: 

1) In the low gain regime {u p l c <C 1), photon pairs are generated uniformly along the crystal. In this case the choice 
Az l c /2 lowers the effect of diffraction and non collincar propagation, since the mean propagation distance that 
photon pairs must undergo to reach the imaging plane tt from the point in which they are created is minimized. The 
factor depending on the refractive indexes takes into account that photons propagate in a dense medium rather than 
in free space. 

2) On the other side, the fields imaged from a generic plane inside the crystal lying at distance Az from the output 
face can be obtained by a virtual free space back-propagation from plane z — l c to plane z — l c — Az. As already 
mentioned, free space propagation leads to an angular divergence of the signal and idler beams with a mean aperture 
angle equal to etc ~ P2 along the walk-off direction. The signal and idler photons are therefore pull apart a distance 
Ay = olqAz w p 2 Az when observed in plane z — l c — Az. In order to compensate this effect the two detection 
areas must be separated by the same distance along the the y-axis. This explain why the value of Ay for which the 
fluctuations are minimized in a given imaging plane is proportional to Az, as it appears from the cigar shaped region 
of low fluctuations along the diagonal direction illustrated in the contour plot of Fig. 10. In particular for Az = Az opt 
walk-off is compensated by taking Ay — acAz opt = pil c l2 and we obtain thereby the shifts given by Eqs. (48a) and 
(48b) for the limit a p l c < 1. 

3) The factor depending on the gain parameter in Eq.(48a), which decreases as a p l c increases, can be understood 
by noting that in a high gain regime most of the photon pairs are generated in the last part of the crystal, because 
of a cascading effect. Hence in order to minimize the propagation distance from the point where they are created to 
the plane n, this plane should be taken closer and closer to the crystal exit face as the gain is increased. 

Thank to this procedure, the noise in is considerably lowered with respect to measurements performed at plane 
z = l c with aligned detection areas (i.e. with Az = Ay = 0). Although the values (48a) and(48b) have been evaluated 
within the plane-wave pump approximation, the numerical simulations demonstrate that the procedure works well 
even when the pump has a finite size. It should also be emphasized that both shifts (48a) and (48b) proved to be 
equally necessary in order to obtain this improvement. 
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VII. CONCLUSIONS 



The results of this paper demonstrate that SPDC is able to display spatial correlation effects at the level of 
quantum fluctuations even in a regime of high gain, i.e. when the down-converted photons form macroscopic fields. 
The quantum origin of the S/I correlations, which lie in the position and momentum entanglement of the twin photons 
building up the two beams, can be best demonstrated in a type II phase-matching configuration where both near field 
and far field measurements can be implemented. We showed numerically that spatial far field correlations of quantum 
origin are observable when the pump beam waist is in the millimeter range and the detection areas are larger than the 
resolution area of the system. Near field correlations seem more difficult to observe experimentally since propagation 
tends to destroy the position entanglement of the generated photon pairs. We proposed a detection scheme which 
allows to optimize their measurement by compensating the detrimental effect of diffraction. These results are strongly 
related to a recent paper of ours [27], which discusses the topic of entangled imaging and extends this technique to 
the macroscopic domain. The simultaneous presence of spatial entanglement in both the near and the far field plays 
a crucial role in the analysis of [27]. In this paper we provide a more quantitative analysis of the level of quantum 
correlation which is present in the far field on the one hand and in the near field on the other. 
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APPENDIX A: INPUT-OUTPUT FORMALISM 



The finite pump pulse bandwidth in space and in time generates coupling between all modes of the S/I field, 
deteriorating thereby the perfect correlation between the (q, 0)i and (— q, — Q,) 2 modes. In this section we generalize 
the input-output formalism of the PWPA to the finite pump case. Input-output transformations (14) are replaced by 
the more general linear transformation 



ai(z,q,Ct) 



dq' I dft' U^z-q^-q Tl')a\ n (q + q ',0 + 0') + Vi(z;g,0;q ',9,')^{-q + q ',-0 + 0') , (Ala) 



a 2 (z,q,Q) = J dq ' J dn' [v 2 (z;q,Q;q ', 0'K nt (-9 + <7 ',-0 + 0') + U 2 (z;q,Q;q '0> 2 "(q + <f ', O + O')] , (Alb) 

which express the fields in a generic plane z inside the crystal in the form of a convolution integral with the input field 
Fourier modes a*"(<f, O) = a,j(z — 0, q 7 O), j — 1, 2. ^From Eqs.(7), we can obtain a fully equivalent set of equations 
for the propagation kernels 



dU!{z;q,n;q ',0') 



dV 2 (z;g,0,g',0') 
8z 



idi (q, VL)Ui (z; q, O; q" , O') 
dO" 



(A2a) 



A (z, q", 0")V 2 *(z; q " - q, O" - O; q ' + q ", O' + O") , 



iS 2 (q,n)V 2 (z;q,n,q',n') 



(A2b) 



+ae' 



-iAoz 



dq" 
~2V 



dfl" 



A (z,q", n")U* (z; q " - q, O" - O; q' + q ", O' + O") . 



The equations for the remaining kernels Vi and U 2 can be obtained by interchanging indexes 1 and 2, and the following 
initial conditions must be fulfilled 



U j (z = 0;q,n;q',n') = S(q')6(n'), 
Vj(z = 0;q,n;q ',0') = , 



(.7 = 1,2). 



(A3a) 
(A3b) 



It can be shown that the solutions of this set of equations satisfy the relations 
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J dq" j dSl" [Ui {z;q,Q;q" - q, Q" - fi)W* (z;q,Q;q" — q',Sl" — SI') (A4a) 

-Vi(z; ft ft; g " + £ ft" + njVf^; q, SI; q" + q', Si" + SI')} = W - <f)K^ - «') , 
y d?"y dO / 'Wi(^;g,0;g"-g,n"-n)V 2 (^;g,0;g" + g",0" + n') (A4b) 

= y rf( ?"y df2''Vi(2;g,n;9'' + g,0'' + n)W 2 (-2;g,0;g''-g',0''-0') , 

which generalize the unitarity conditions (17) beyond the case of plane-wave and cw pump. 

The normally ordered photon number correlation function in a generic transverse plane z can be written as 

G^\x,t, x ' ,t') = (at (2, x, t)ajj(z,x ' , t')a,j(z, x ',t')a,i(z, x, t)) - (a\(z,x,t)a,i(z, x, t))(aj(z, x ',t')a,j(z,x ' ', t')) (A5a) 

= \(4{z,x,t) aj (z,x' ,t'))\ 2 + \(a i {z,x,t)a j (z,x' ,t'))\ 2 . (A5b) 

In the last identity we made use of the general property characterizing fields with a gaussian statistics, which allows to 
write the fourth order field correlations as a sum of products of the second order correlation functions (see e.g. [22]). 
We now consider explicitly the transformations relating the S/I fields in the planes where detection is performed, to 



those on the crystal output face, a° ut {q,t): 



aj {z, x,t) = J dq hj(x, q)a° ut (q, t) = J dx ' hj(x, x ')a° ut {x ', t) (j = 1,2) 



with 



dr 1 , 
hMQ)= I ^ q ' x h^x') 



(A6) 



(A7) 



Propagation outside the crystal can include several optical devices, such as lenses and polarizing beam splitters, and 
can take different paths for the signal and idler beams (see scheme of Fig. 4). We shall assume however that it occurs 
without losses and this latter condition implies that the Fresnel kernels hj (x, q) satisfy the relation 



dx h*(x,q)hj(x,q') = 6(q- q) , (j = 1,2) 



(A8) 



which can be obtained by requiring that the commutation rules (4) are preserved in the transformation. When 
the detection time Td is much longer than the coherence time Qq , the photon number self- and cross-correlations 
measured over two detection areas i?i and R 2 , as defined by Eq.(30), can be written as 



dSl [ dSl' [ dx [ dx' f dq f dq'h^x^h^x ' \q'){a^ n (q,Sl)l 



(5N 1 5N 2 ) Z = dSl dSl' 



dx I dx ' 

Ri 



dq / dq'h^X, q)h 2 {x ', q ')(€.?* {q, n)a%* (q ', ft')} 



(A9a) 
(A9b) 



where we used relation (A6) in order to express the second order field correlations appearing in Eqs.(A5b) in terms 
of the output field operators a° ut (q,Sl), j = 1,2. The correlation functions of the output fields can be expressed in 



terms of the propagation kernels defined by Eq.(Ala) evaluated at plane z = l c as 

(a? u * t ($;ft)a; ut ($",n / )) = J dx j die *(?-? ')-*-i(n-n)* y*(/ c ; «f, fi; x, i)Vi(/ c ; <f \Sl';x,t) 
(af^ftfija^^'.fi')) = J dx J dte i«4-?')-*-*(n+n)* - Q . - t )v 2 (/ c ; q', ft'; x, t) 



with 



Uj(z;q,Sl;x,t) = J ^ J 



dt 



■x—ifl r t 



2ir 



Uj(z;?,Sl;?',tf) (.7 = 1,2). 



(AlOa) 

(AlOb) 
(AlOc) 

(All) 



A similar definition holds for the functions Vj(z; q, ft; x, t), j = 1, 2. In case R\ and R 2 intercept all the photons of the 
signal and the idler fields generated in the down-conversion process, it can be shown that the variance of N- = N\ — N 2 
vanishes if condition (A8) is fulfilled (i.e if free propagation occurs without losses). However, in general such a result 
does not hold if the two detectors collect photons only from finite portions of the two beams. 
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APPENDIX B: APPROXIMATE SOLUTION IN THE QUASI-STATIONARY REGIME 

In this appendix we derive an approximate analytical solution of the propagation equations (7) assuming the spatial 
and temporal frequency bandwidths of the pump are small but finite; more precisely we assume that conditions (24) 
are satisfied. 

It is useful to write the propagation equations in the reference frame comoving with the pump field envelope, 
whose coordinates are related to the original laboratory coordinates through the linear transformation t' = t — k' z, 
y' = y + pqz. In Fourier space this correspond to multiplying the Fourier components of the S/I field envelopes by 
e -i[k n+p q y ]z ^ jyf orc p rec i S ely we consider the transformation 

afaftil) = A^-W-wKjiz,?,!!) , (j = 1,2) , (Bl) 

where the constant phase factor e l ~ z has been added in order to eliminate e' 1 aZ from the convolution term in 
Eqs.(A2). Input-output transformations (Ala) and unitarity conditions (A4) still hold for the transformed kernels 

U/(z;q,n;q\n , )=e^- k '» n - poq ^ z U 1 {z;q,n;q',n') , (B2a) 
V 3 '(z;q,n;q-i,n>) = e i \Q~ k '° sl ~ Poqy \ !l Vj(z; q, SI; q ', SI') , (j = 1,2), (B2b) 
which satisfy the propagation equations 

(B3a) 

0, q ", Sl")V' 2 *(z; q" - q, Si" - 0; q' + q ", Si' + Si") , 

(B3b) 

0, q ", Sl")U[*{z; q " - q, O" - O; q' + q ", O' + Q") . 

The explicit form for the propagation of the pump field (9) has been used and we defined the new detuning parameters 

6' 3 (q, 0) = ^ + Sj (q, D.) - p q y - k' Q (j = 1, 2) , (B4a) 
S' (q,n) = \k>^ ^-f . (B4b) 

Remembering the hypothesis that the pump envelope Fourier transform at plane z — has the gaussian form (11) 
with Sqo <C qo and 6ojq ri , we can now apply the following approximations: 

1) The phase term in the convolution integrals can be neglected, since 5' (q, Cl)l c is at most on the order of 
max |(5g§/(jQ, SuJq/£Iq} in the region where the Fourier transform of the pump envelope is not negligible. 

2) We expect that the propagation kernels defined by Eq.(Ala) are characterized by the slow variation scale (go, ^o) 
in their unprimed arguments, while they are strongly peaked in the origin of the primed variable space (q ',&'), in 
which they have a much faster variation scale. This assumption is justified by the form of the solutions of Eqs.(B3) 
obtained in the PWPA limit. As Sqo/qo — > and £wo/^o — * 0, the primed propagation kernels satisfying initial 
conditions (A3) take indeed the simple form: 

U' 3 (z;q,n-q\n r ) = 5{q ')5{Q.')A^- k '^- poq y\ z U j {z; q, SI) , (B5a) 
V 3 (z;q,n-q ',fi') = 5{q r )5(n')e l i^ l - k '° n - poq ^ z V : j(z; q,fl) (B5b) 

where the gain function Uj(z; q, SI) and Vj(z; q, ft) are given by Eqs.(15), with l c being replaced by the ^-coordinate. 
This latter hypothesis allows us to neglect the dependence on q " and SI" in the first argument of the kernels in the 
convolution integrals at the r.h.s. of Eqs.(B3a) and(B3b). With these approximations, we obtain a system which can 
be solved analytically and which acquires its simplest form when written for the kernels Fourier transformed in their 
primed arguments (see definition (All): 



9Wi'(z;q,0;q \Sl') 
dz 



dV 2 , (z;q,n,q',n') 

dz 



iS'^Slpx' \z;q,Sl-q',Sl') 

J 27T J V27T 

i6' 2 (q,n)V 2 '(z;q,Sl,q',tt) 

+ [^[^e<^ z Mz 
J 2tt J V2^ 
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<9ZV(<f,ft;x, t) 

dz 

dV 2 '(q,Sl,x, t) 

dz 



= i8[{q,Sl)Ux'{q,Sl;x, t) 

+aA Q (z = 0,x,t)V 2 *(-q, -Sl;-x, -t) , 
= iS' 2 (q,n)V 2 '(q,n,x,t) 

+aA (z = 0,x,t)U^*(-q,-Sl;-x,-t) , 



(B6a) 
(B6b) 



The solution of this system satisfying initial conditions (A3), which now read U j(z — 0;q,Sl;x,t) = l/(2n) 3 / 2 , 
V'Az = 0; q, SI; x, t) = 0, [j = 1, 2), are 



U\ '(z; q, SI; x, t) = exp 
Vi '(z; q, SI; x, t) = exp 
U 2 '{z; q, SI; x, t) — exp 
V 2 '(z; q, SI; x, t) = exp 



. S[(q,Sl)-5' 2 (-q,-Sl) ' 
i x z 





2 




S[(q,Sl) 


-S' 2 {-q, 


-S}) 




2 




S' 2 (q,Sl) 


~S[(-q, 


-SI) 




2 




S' 2 (q,Sl) 


-5[(-q, 


-SI) 



U{z;q,Sl;x,t) , 
V(z;q,Sl;x,t) , 
U (z; —q, —SI; —x, —t) 
V(z; —q, —SI; —x, —t) 



(B7) 



with 



U(z; q, SI; x, t) = 
V(z; q, SI; x, t) — 



(2tt) 3 / 2 
1 



coshr(<f, SI, x, t)z + i — y=rj!^l sinhr(<f, SI, x, t)z 



aA (x, t) 



(2tt)3/2 T(q,Sl,x,t) 



2T{q,Sl,x 1 1) 
sinhr(g, SI, x, t)z , 



F(q, 0; , x , t) = y o*A*{2, t) , 

A(q, SI) = 6[(q, SI) + 5' 2 (-q, -SI) = A + S^q, SI) + 5 2 {-q, -SI) . 



(B8) 



(B9) 
(BIO) 



Clearly, as 5qo,Suj — > these functions function looselose their dependence on the space-time coordinates (x, t) and 
we recover the PWPA solution expressed by Eqs.(15) and (16). In the more general case in which the ratio Sqo/qo 
and Sqo/qo are small but finite, the self-correlation function (AlOa) is peaked at (? = q, while the cross-correlation 
function (AlOb) is peaked at <f = —q, the width of both peaks being on the order of Sqo- This is more clearly seen 
by considering the special limit in which l c — > and the gain parameter aA p l c remain a finite quantity. In this case 
the effects of linear propagation, such as diffraction, dispersion and walk-off become negligible and the propagation 
kernels (B8) lose their dependence on q and SI (as it can be inferred from the fact that the characteristic bandwidths 
defined in Eq.(20) go to infinity as l c — > 0). The correlation functions (AlOa) and (AlOb) are then simply the Fourier 
transforms of sine and cosine hyperbolic function of aA (x, t)l c , calculated in q — cf and in q + <f respectively. We 
also verified that that these approximate solutions looks very similar in shape to those obtained with the complete 
numerical model (see Fig. 6), although for the chosen crystal length of 4 mm propagation effects are far from being 
negligible. 

Considering the specific case of a far field measurement with the / — / lens system described in Sec.V, the Fresnel 
kernels defined by Eq.(A6) take the form hi(x, q) = h 2 (x,q) = —jj-5 (^q— jjX^J , as it can be verified by substituting 
expression (21b) into Eq.(A7). Eqs.(A9) reduce then to 



(: (dm) 2 :)„, = (: («5iV 2 ) 2 :) n > = dSl dSl' dq / dq ' 

J J JQi JQi 

(SN 1 5N 2 ) 7T/ = [ dSl [ dSl' [ dq [ dq' [ dq [ dq , (a° 1 ut (q,Sl)a 2 mt {q",Sl' 
J J JQx JQ 2 J J 



dq\ dq (aT\q,Sl)aT\q',Sl')) 

2 



)) 



(Blla) 
(Bllb) 



where Qi and Q 2 indicate the regions in the spatial frequency plane corresponding to the two symmetrical detection 
areas Ri and R 2 according to the mapping x —> jjx. Knowing that ((5N-) 2 ) —> as R\, R 2 —> oo, the localization 

of the cross- and self- correlation functions (AlOa) and (AlOb) on a area on the order of Sq^, for <f = q and <f = —q 
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respectively, guarantees that nearly complete noise reduction is achieved if the area of the two detectors is large 
compared to the resolution area determined by the pump beam waist, that is Sdiff — (\f /2ir) 2 5qQ. 
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FIG. 1. Scheme for the observation of spontaneous down-conversion in the far field zone (a). The lens (not shown in the 
figure) is located at z = l c + /. (b) and (c) display the phase-matching curves (25) in the spatial frequency plane for a type II 
(b) and a type I (c) crystal respectively. The symmetrical black square Ri and R2 indicate the location of the detectors from 
which maximal signal-idler correlation can be measured. 





Ao>0 




FIG. 2. Typical far field pattern from spontaneous down conversion in a type II crystal, assuming observation is performed at 
the degenerate frequency (i.e. at Ai = A 2 ). They are obtained for decreasing values of the collinear phase mismatch parameter 
Ao, which makes the radius of the rings shrinks to zero. 




(b) 



A >0 



FIG. 3. Far field pattern in a type I crystal at the degenerate frequency for collinear (a) and non-collinear (b) phase-matching. 
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FIG. 4. Detection scheme to measure spatial correlations in the near field. A polarizing beam splitter (PBS) separates the 
S/I beams. Their near fields, at the plane it : z = l c — Az, are imaged by two lenses (L and L') onto the pixel detectors Ri and 
i?2, which lie in the plane conjugate to plane it. Az and Ay indicate the spatial shifts applied to the optical devices that are 
necessary to to optimize the measurement. 



NEAR FIELD 



FAR FIELD 




FIG. 5. Near field (left) and far field (right) patterns obtained for 5qo/qo = 0.05 (corresponding to a pump beam waist 
if o = 920/im in the LBO case) (a), Sqo/qo =0.1 (wo = 460/im) (b), and 5qo/qo = 0.3 (wo = 150/im) (c). AoZ c = 13.6 and 
<j v lc — 3.0. 
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FIG. 6. Far field correlations: {: 8N(x)SN(x ') i}^/ is plotted as a function of x, for increasing values of the ratio Sqo/qo- In 
(a) the mean photon number distribution profile is also shown (dotted line), x' = 3.6 xo is kept fixed in the region of maximum 
gain. The other parameters are the same as in Fig. 5. 
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FIG. 7. Far field correlation in type I: ((<5A r _) 2 ) 7r //(A r + )^/ is plotted as a function of the detector size for different values of 
8qo/<lo- The other parameters are the same as in the previous figure. 
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FIG. 8. Far field correlation: the ratio {(5N-) 2 )^i / (N + ) w i is plotted as a function of the detector size d for increasing value of 
the ratio Sqo/qo- The parametric gain is a v l c — 4. The negative value of the collinear phase mismatch, Ao/ c = —qc/qo = —74.4, 
is such that the radii of the rings qa vanishes. 
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FIG. 9. Near field correlation: the ratio {(SN-) 2 )^ / {N + 
the pulsed gaussian pump are wo = 332^im (5qo/qo = 0.1) 



is plotted as a function of the detector size. The parameters of 
and to = 1.5 ps (Sujq/Qq = 1.14); the gain is a p l c — 3 and qn = 0. 
The simulations performed applying diffraction and walk-off compensation (squares) are well below the one performed without 
optimization (white triangle). The dashed line corresponds to the analytical solution obtained in the PWPA, given by Eqs.(29) 
and (41) 
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FIG. 10. Near field correlation: contour plot of the ratio {(5N-) 2 }„/{N+} Tr in the (A«, Ay) plane calculated according to 
the PWPA theory. The minimum in the center corresponds to the values given by Eqs.(48). 
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